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We present a quantum algorithm to prepare tire thermal Gibbs state of interacting quantum 
systems. This algorithm sets a universal upper bound D" on the thermalization time of a quantum 
system, where D is the system's Hilbert space dimension and q < | is proportional to the Helmholtz 
free energy density of the system. We also derive an algorithm to evaluate the partition function 
of a quantum system in a time proportional to the system's thermalization time and inversely 
proportional to the targeted accuracy squared. 

PACS numbers: 



The partition function at inverse temperature /? is of 
central interest in statistical physics because all thermo- 
dynamic variables such as energy, heat capacity, entropy, 
and correlation functions can be derived from it. For in- 
teracting systems, the partition function involves a sum 
over an exponential number of terms, so it cannot be 
computed exactly except in few special cases. Thus, ap- 
proximation schemes are required to estimate these phys- 
ical quantities numerically. 

Monte Carlo is often the method of choice to evaluate 
thermodynamic quantities. When the system is classical, 
system configurations a; can be generated randomly ac- 
cording to the thermal distribution Pr(a;) = ^ z{f3)'' by 
a suitable Markov chain. Sampling from this distribu- 
tion provides an estimate of the partition function and 
other quantities of interest. The running time of the al- 
gorithm is determined by the rate of convergence of the 
Markov chain and the number of samples needed to lower 
the statistical fluctuations below a desired accuracy. The 
method of ref . [17] makes it possible to prepare a purifi- 
cation of the thermal distribution in time that scales like 
the inverse of the square root of the spectral gap. Build- 
ing on this work, one of us has recently shown that a 
quantum computer also achieves a quadratic speed-up 
with respect to relative accuracy for the task of evaluat- 
ing partition functions [19 . 

Monte Carlo can also be employed for quantum sys- 
tems using standard mapping between iV-dimensional 
quantum systems to TV -I- 1-dimensional classical systems 
(the N + 1th dimension corresponding to inverse temper- 
ature). However, frustrated and fermionic systems are 
affected by the sign problem jSj which leads to an ex- 
ponential blow-up of the statistical fluctuations of the 
sample. For this reason, Monte Carlo simulations are 
unreliable for these systems. 

This Letter presents two complementary results. First, 
we demonstrate a quantum algorithm to prepare the 
thermal Gibbs state p{P) = ^gjjs) '-'^ locally inter- 
acting quantum system, where Z{P) = Tr{e^'^^) is the 
partition function. The problem of thermalizing a quan- 
tum system with a quantum computer was flrst raised in 



[TH] where heuristic methods were proposed, but no rigor- 
ous upper bounds have been derived (for classical systems 
the problem was flrst raised in [H]). Our results set a 
universal upper bound D" for the thermalization time of 
a quantum system, where D is the system's Hilbert space 
dimension and a < | is proportional to the Helmholtz 
free energy density of the system. 

Second, we present a quantum algorithm to evaluate 
the partition function of quantum systems within relative 
accuracy e. The complexity of the algorithm scales like 
the inverse of (e when the system is classical) and 
proportionally to the thermalization time of the system. 
This algorithm is not affected by the sign problem and 
works for any locally interacting quantum system. It is 
a full quantum generalization of the algorithm presented 
in [T71 [H] for the evaluation of partitions functions of 
classical systems. The combination of our two algorithms 
yields a universal upper bound for the time needed to 
evaluate a partition function with a quantum computer. 

We consider a system composed of n interacting d- 
level particles with a k- local Hamiltonian H — hj 
where each term hj has bounded norm and acts on at 
most k particles (fc is a constant). The most physically 
relevant case is of course k — 2. It is one of quantum 
information science's most celebrated result that /c-local 
Hamiltonians can be simulated efficiently on a quantum 
computer. More precisely, for those Hamiltonians, the 
time evolution operator U{t) = e~*^* can be approxi- 
mately expressed as the product of a sequence of dis- 
crete one- and two-qubit gates, and the number of gates 
scales essentially linearly with the duration t of the sim- 
ulated process and polynomially with the system size n 

Given the system's Hamiltonian and a temperature 
1//3, we call a thermalization process any (possibly time- 
dependent) local Hamiltonian on the system and a bath 
(of size scaling polynomially with the number of particles 
in the system) that generates, after time T, a unitary 
transformation V such that 

Trbath \v |0)(0|bath) p(/3), (1) 
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or more generally that this equation holds within some 
accuracy 6. The system's thermalization time rth(/3) is 
the duration of the shortest thermalization process. Note 
that in our definition, the temperature dependence is en- 
coded in the details of V instead of in the initial state 
of the bath as one would expect physically. This is to 
rule out trivial thermalization processes where the bath 
is prepared in the desired state and V simply swaps 
the system and the bath. Nevertheless, the temperature 
dependence of the thermalization process we will present 
could equally well be encoded in the initial state of the 
bath with a simple modification to our algorithm. 

Thermalization processes that occur in nature consti- 
tute of local interactions (all forces of nature are two- 
body) , so can be efficiently simulated by a quantum com- 
puter, i.e., in a time proportional to the duration of the 
thermalization process [TJ [3l US [2^ . In what follows, we 
present a universal thermalization process and character- 
ize its complexity. Then, we show how to use this~or any 
other-thermalization process to evaluate partition func- 
tions with a quantum computer. 

For sake of analysis, we shift the system's Hamilto- 
nian by a constant E* to ensure that H is positive. 
This changes the partition function by a factor e~^^ 
that we keep track of implicitly: we henceforth assume 
< H < E'max where i?max is a known upper bound of 
the Hamiltonian (since H is fc-local, i?max scales polyno- 
mially with the number of particles and can be straight- 
forwardly upper bounded by triangle inequality) . Wc de- 
note the eigenvalues and eigenvectors of the Hamiltonian 
H\a)^Ea\a). 

We will make extensive use of the system's time evolu- 
tion operator with t — 7r/(4ii^max), that we simply denote 
U. As mentioned above, U cannot be simulated exactly 
with a quantum computer. Instead, we can produce U 
that is a good approximation to it. This will unavoidably 
limit the accuracy of our thermalization process and of 
the estimated partition function. Indeed, defining the 
effective Hamiltonian H — —jLogU (by definition of t, 
there is no multi- value problem) , our algorithms will pre- 
pare the Gibbs state p{(3) and evaluate the partition func- 
tion Z{P) associated with H, not H. However, we show 
in Appendix A that an accurate simulation ||?7 — [/|| < e 
leads to an accurate effective Hamiltonian \\H — H\\ < 
KE^ii^e (where is a small constant). This, in turn, 
implies an accurate estimate of the partition function 
\Z{(3) - i(/3)| < 2K(3Erna^eZ{(3) (Appendix B ) and a 
high fidehty Gibbs state F[p{P),p{P)] > 1 - K' (}E,^^^e 
(Appendix C ). These bounds can be understood intu- 
itively by taking a first order expansion in the small (3 
limit. 

Because the exact time required to simulate U within 
accuracy e depends on the details of the system, we mea- 
sure the complexity of our algorithms in terms of the 
number of times they need to implement U . Hence, the 
true running time of our algorithms have an additional 
arbitrarily weak dependence [3^ on the targeted accuracy 
e of the partition function and Gibbs state. 



Our algorithm requires a cooling schedule = /3o < 
f3i < ... < Pi ^ (3 such that Fk := ^zi}T)^ ^ 1/2 for 



all k. Writing Fk - -^jj^ 
where A/3j. = /3fc+i — we see that a polynomial- length 
cooling schedule can always be constructed by choosing 
A/3fc = log(2)/i?i„ax, which implies a cooling schedule of 
length I = ii^niax/?/ log 2, but shorter schedules can be 
used in most cases. Notice how the choice of a lower 
bound on the Hamiltonian E* affects these ratios. For 
this reason, it is desirable that —E* be as close as possible 
to the true ground state energy of the system to decrease 
the length of the cooling schedule. 

The partition function is expressed as a product 



Z(/3fc) = Z(/3o)FoFi 



■ F^ 



(2) 



and will be computed by evaluating each fraction and 
using the fact that -Z(/3o) = the dimension of the sys- 
tem's Hilbert space. If each ratio Fj. is evaluated within 
accuracy |, the resulting relative error on -Z(/3) will be 
0(t). Our universal thermalization process requires eval- 
uating the Fk sequentially: the fraction F^ is evaluated 
from the Gibbs state p{Pk), and the value of Z{f3j.) is 
needed to prepare p(l3k+i). 

Our method makes use of quantum phase estimation 
(QPE) [10l[T5]. This is a transformation that operates on 
the "system" register and an m-qubit "energy" register 
initiahzed to the state |0m) as follows: 



10. 
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Here, W denotes the Hadamard transform VF|Om) = 
2~^ Y^r^o^ k) and FT denotes the Fourier transform. 
The central gate applies r repetitions of the time evolu- 
tion operator U to the system register, where r is dic- 
tated by the state of the lower register. Since r can be 
as large as 2™ — 1, the running time of this procedure 
is 2"*. The label E runs over the discretization of the 
interval [0, 8E^^^), i.e. E e x {0, 1, . . . , 2™ - 1}. 

The function f{E, Eg) is the discrete Fourier transform 
'', so it is highly peaked around the value 



of g{r) : 
E = Ea, i.e. 



itE„ 



\fiE,Ea) 



sin(^f2") 



Thus, QPE essentially measures the energy of the sys- 
tem and writes it down on the energy register, up to 
some fluctuations associated to the width of the function 

f{E,Ea). 

Let us first describe our algorithm to prepare Gibbs 
states by ignoring these fluctuations and return to them 
later. In fact, our algorithm will prepare purified marked 
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Gibbs states 



(3) 



that are defined on one system register, one scratchpad 
register, and one energy register. The term "purified" 
refers to the fact that the state is pure and has a 
reduced system density matrix equal to the Gibbs state 
at inverse temperature Pk- The term "marked" refers to 
the fact that the energy register encodes the energy of 
each system eigenstate. 

To obtain the state |/3o), we must first prepare the sys- 
tem and scratchpad in any maximally entangled state 
-D^2 |a) (g) and then apply QPE to the system 
register (ignoring fluctuations of /). The finite tempera- 
ture quantum Gibbs state \(3k) is obtained from |/3o) and 
an ancillary qubit initially in the state |0) that we rotated 

by an angle 9{E, Pk) = arcsin(e' 
energy register, resulting in 



2 ) conditioned on the 



\^k) = D---Y,\a)®\4>a)<i^\Ea)(E\e{Ea,f3k)) (4) 
l/3fc>|0) + ... (5) 



|0). Indeed, applying q = 
of two reflections (/ — 21^*^ 



D 

where \0) = sin^^|0) -|-cos0|l) and the ellipsis represents 
terms in which the ancilla qubit is in the state Note 
that the required conditional rotation can be expressed 
as the product of m rotations of angle arcsin(e t™^) 
controlled by a single qubit. 

We can use Grover's algorithm |6] to amplify the over- 
lap of Y^k) with the projector Hg = J (g) |0)(0| associated 
to the subspace where the ancillary qubit is in the state 

-2p^J times the sequence 

«'fe|)(l - 2no) to the state 
Y^k) will boost the amplitude of the term |/3fe) (g) |0) from 

its initial value \J ^''^''^ to nearly 1. However, q cannot be 
computed exactly because at this stage of the algorithm, 
Z{(ik) is known only within relative accuracy This 
is nonetheless sufficient to amplify the amplitude of the 
term \(3k) ® |0) to 0(1). An amplitude 1 - 0(f) can be 
achieved using fixed-point search [Zl H] at an additional 
multiplicative cost of log ^ . 

The fluctuations of f{E,Ea) will in general invalidate 
the procedure we have described. If we perform QPE 
with a slightly larger number of qubits — i.e., if we use an 
(to -I- 4)-qubit energy register — the probability that the 
estimated energy E deviates from its true value Ea by 
more that 2^'"/t will be less than 1/16 [15 . To further 
suppress these fluctuations, we can perform rj indepen- 
dent QPE procedures and compute the median AI of the 
77 results. The probability that this median value devi- 
ates from the true energy by more than 2~™/< is less 
than 2"'' [M]. 

The evaluation of the median can be done by a coher- 
ent quantum process. In other words, there is an efficient 



quantum circuit on rj energy registers and one TO-qubit 
"median" register that maps . .(g)|i?(''))g) 

|0™) to 1^(1)) (g|i;(2))(g...(g|s(''))(g|M(^(i), £;(!),...)). 
From the result stated above, we know that on input 
eigenstate \a), only two median register states \E^) and 
\E^) encoding the closest TO-bit estimates of Ea can have 
an amplitude of magnitude larger than \Jl~'^ . This can 
be used to prepare a "good enough" version of the infinite 
temperature purified marked Gibbs state 



(6) 

where [a^p + |/3,|2 > l - 2-", |S± - Ea\ < 2-^/t, 
and I bad) represent the components where the median 
deviates from the true energy by more than 2~™/i, so 
(badjbad) < 2^''. The scratchpad states contain 
the original states \4>a) used to purify the infinite temper- 
ature Gibbs state and the rj (possibly entangled) energy 
registers used for the computation of the median, whose 
value is encoded in the third register. 

With TO = log2{f3i/te) the error caused on the Gibbs 
state by the roundoff of the energy estimate will be 0(f ) 
(Appendix C). Similarly, this roundoff error will lead to 
a 0(f) relative error on our estimate of Fk (see Ap- 
pendix B). The error due to the |bad) component of 
Eq. ([g]) requires a bit more attention. The reason is 
that the relative weight of the good and the bad com- 
ponent [i.e., first and second term of Eq. (|6|] is in gen- 
eral not preserved by Grover's amplification. Although 
the norm of [bad) starts out small, it can increase by a 
factor e''''^'"''^/^ during the amplification process. Set- 
ting 77 = [ln(j) + /3fci?,„ax]/ln(2) insures that the norm 
of [bad) will remain bounded by | after the amplifi- 
cation process, so can safely be ignored. This com- 
pletes our presentation of the universal thermalization 
process, and demonstrates that the Gibbs state can be 
prepared with accuracy (5 = | by a local process in time 

e O ^^f^!^ log i[l0g i + /3i?,nax]^ 

We now describe how the state |/3fe) can be used to 
estimate the ratio Fk . Similarly to the procedure that led 
to Eq. ([5]), we append to the quantum Gibbs state \(3k) 
yet another ancillary qubit in the state |0) and rotate 
it by an angle d{E, A/3k+i) conditioned on the energy 
register, resulting in 



E 



ZiPk) 

Fk\f3k+im 



a)(g,\4>a)(E\Ea)(E>\9iEa,APk+i)) 
+ ■•■ (7) 



where the ellipsis represent terms in which the ancilla 
qubit is in the state |1). Clearly, the squared norm of 
Ho|<I>fe) is equal to the quantity Fk we are trying to esti- 
mate. Thus, Fk can be estimated by quantum counting 
[4], i.e. using phase estimation to estimate the eigenvalue 
of the unitary matrix (J — 2|<I>fc)(<i>;;|)(l — 2Ho) associated 
to the Jordan block that supports An accuracy 
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FIG. 1: (Color online) Scaling exponent as a function of the 
inverse temperature for the one dimensional Ising model de- 
fined by the Hamiltonian H{J, g) = JI^^ gcri'' + Jai'Vi*'''^-' for 
different values of g/ J (g/J=l is critical). The absolute val- 
ues of g and J are fixed by \\H{g, J)\\ = \\H(Q, 1)||. Like the 
partition function, a is symmetric under J ^ g. 

I requires j uses of the circuit that prepares so 
the total time required to evaluate the partition func- 
tion with accuracy e with constant success probabihty is 

^ (\/Sy^"^)- 

More generally, any thermalization process can be used 
to evaluate F^.. Indeed, a purified Gibbs state can be 
prepared by simulating the thermalization process with 
a quantum computer, substituting the system's maxi- 
mally mixed state in Eq. ([T]) by a maximally en- 
tangled state X^o I'*) ^ I'f'a)- Hence, this requires a 
time proportional to the thermalization time of the sys- 
tem. This state can be marked using multiple QPEs and 
evaluating their median as described above. An error 
at most e is achieved with constant probability in time 

The complexity of our algorithm scales with the sys- 
tem size as \J — D"' where D — d/^ is the Hilbert 

space dimension and we can express the scaling exponent 
a = ■^h{(3) j log((i) as a function of the Helmholtz free en- 
ergy density /i(/3) (energy is defined relative E*). Figure 
[l] shows this exponent for the Ising model as a function 
of inverse temperature. For temperatures of order unity 
in natural units, this scaling parameter approaches ^, 



which naturally coincides with the running time found 
in [16] to compute the ground state energy of the sys- 
tem. At higher temperatures however, the algorithm can 
be significantly faster. Surprisingly, the running time is 
shortest at the critical field gj J — 1, where the correla- 
tions in the system are the strongest. 

For a classical Hamiltonian, the energy of a given con- 
figuration can be computed exactly. In that case, our 
thermalization process produces a Gibbs state with ac- 
curacy (5 in time \J 2.{j3k) ^'^S | ■ The only approximation 
to the partition function involved in that case comes from 
quantum counting, so Z{I3) can be evaluated with accu- 
racy e in time O {^j. \ z{j3k) ) ' ^^^^ should be com- 
pared to the time O ^^j. \ ^j=^ required by the quan- 
tum simulated annealing algorithm |171 119] , where 5k is 
the spectral gap of a Markov chain whose fixed point is 
the Gibbs state of the system at inverse temperature (3k ■ 
Thus, our method offers a speed-up for the computation 
of a classical partition function whenever no such Markov 
chain can be found with 5 > Z{(3)/D. 

Conclusion — We have presented a quantum algorithm to 
prepare the Gibbs state of a quantum system in a time 
that grows like the ath power of the system's Hilbert 
space dimension, where a < 5 is proportional to the 
Helmholtz free energy density h{f3) at inverse tempera- 
ture p. This sets a universal upper bound on the thermal- 
ization time of quantum systems. However, our universal 
thermalization process fails to recognize special proper- 
ties of the system-e.g. large energy gap, absence of long 
range correlations-that could potentially speed up the 
computation. We are actively investigating this problem. 

We have demonstrated how the ability to thermalize a 
quantum system on a quantum computer leads to an al- 
gorithm to estimate its partition function with accuracy 
e in a time proportional to the system's thermalization 
time and e~^. Our method provides a complete quantum 
generalization of the work reported in |17l 119] for parti- 
tions functions of classical systems and may provide an 
additional speed-up when there exist no rapidly mixing 
Markov chain to prepare the system's Gibbs state. 
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APPENDIX A: SIMULATION ERROR FOR THE 
EFFECTIVE HAMILTONIAN 

Let H he a, self-adjoint operator with \\H\\ < 7r/4t 
and let U = e**-^. Denote by \a) and Xa the eigenvalues 
and eigenvectors of H, respectively. Using techniques for 
simulating Hamiltonian time evolutions, we can realize 
a unitary U with \\U — U\\ < e. We prove that there is 
a self-adjoint operator H such that \\H — H\\ < ^ and 
^itH _ jj TffiyQYe K is some small constant. 

Let := C \ (— oo, 0]. We use the logarithm 

logz = log jzj -t- i argz , (Al) 

on Ctt , where arg z is defined to be the unique number 
6 e (— TT, vr) with z = |z|e'^. Let 7 be the contour defined 
in Fig. [2] and denote its interior 1(7). 

The spectra of U and U are contained in 1(7). This 
is clear for U since \\tH\\ < 7r/4. Assume to the con- 
trary that U has an eigenvalue e^^ with |A| > 7r/3. 
Let lip) denote the corresponding eigenvector of U and 
Pa = For e < 1/4, this leads to the following 




FIG. 2: For _R > 1, 7 is the contour defined as the join of 
the outer arc 7i(s) = Re'" for s G [n/2,n/2], the upper line 
segment 72(3) — (R — s)i for s € [0,R— 1/R], the inner arc 
73(s) = {1/R)e-" for s G [-7r/2, 7r/2], and the lower line 
segment 74(5) — (—1/R — s)i for s G [0,R— 1/-R]. We denote 
by 7*1 7ij 72 ) 73) and 74 the corresponding images. The 
contour 7 is traversed in the positive sense as indicated by 
the arrows. 

contradiction 

e > MU\yj) - {^PpM (A2) 
= lE^'4e^'"-e^')l (A3) 

a 

> |e*f _ e^f I = 2sin;^ « 0.26 > ^ . (A4) 
- I I 24 4 ^ 

Generalizing this argument and using the inequality 
< 2sin0 for (f> G [0,7r/2], we arrive at the inequal- 
ity \\H -I- cjj < \\H -\- c\\ -I- e/t for all real c. This implies 
that the spectrum cr{H) of H is contained in the interval 
(min{CT(iJ)} - f,max{cr(i7)} + f). This result will be 
used in Appendix [C] 

For our purposes, it is important that the logarithm 
function defined above is Lipschitz continuous in 1(7). 
More precisely, we have 

|loga-log6| < Kja- &| . (A5) 

This is seen as follows. We may assume w.l.o.g. that 
a = \a\ and b = |6|e*^ with 9 e [0, tt]. We have 

I log a - log 6| 2 = (log \a\-log\b\f + 0^ (A6) 

and 

\a-b\^^{\a\-\b\r + \a\\b\4sm^^-. (A7) 

Using the fact that 1/R < \a\, \b\ < R in 1(7), we clearly 
have 

|loga-log6|2 <i?2(|a|-|6|)2 (A8) 

since the real logarithm is Lipschitz continuous with con- 
stant R on the interval [1/-R, 00), and 

e'<^\a\\b\Asm'^^ (A9) 

following from the inequality 4> < 2sin0 on [0, ^]. Thus, 
the Lipschitz constant is k = max{-^, R}. 
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Let A be an arbitrary diagonalizable matrix and S an 
invertible matrix such that 



A = S diag(ai, . . . ,a„) S" ^ 



(AlO) 



Assume that ai, a2, . . . , a„ are in C^r- Then, it is natural 
to define the logarithm of A as follows 

LogA = S'diag(logai,loga2, . . . ,loga„) S*"^ . (All) 

It is possible to extend the definition of the logarithm 
to all matrices A G M„(C) with spectrum <j{A) C 
with the help of the Jordan canonical form of A. This 
gives rise to the so-called primary matrix function LogA 
associated to the scalar-valued stem function logz. We 
refer the reader to [5J Chapter 6 "Matrix Functions" ] for 
more details. 

Since U and U are unitary, and thus diagonalizable, 
we can define their logarithms using the formula in 
eq. (All I. It is clear that 



Logf/ = itH . 

We would like to be able to conclude that ||LogC/— Logf/| 
is small since \\U — U\\ < e is small. This is true when U 
and U commute. In this case, we may assume that they 



are both diagonal. Hence, everything reduces to eq. (A5) 



since it suffices to compare corresponding diagonal en- 
tries. 

To prove the desired result in the non-commutative 
case, we need to consider the following result which is 
established in Theorem 6.4.20 in [9 . Let V = {A e 
Mn{C) : a{A) C 1(7)}. For any A eV, define 



LogA = 



1 

27ri 



{log z){zl - A)-^dz. 
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Then, Log is continuous on V, and 



(A12) 



(A13) 



Let A and B be two arbitrary unitary matrices with 
a{A),a{B) C 1(7). We now show that 



\LogA - Log^ll < K\\A - B\ 



(A14) 



holds for all such unitaries A and B, where K is some 
constant. We have 



I LogA - LogS 
„ 1 



(A15) 



(logz)[(z/ - A)-^ - {zl - B)-^]dz\\ (A16) 



— i log zlzl - Ay [A - B){zl - By ^ dz 
«27r 



(A17) 



< 



^ £ I logzl \\{zl - A)-i|| ||A - B\\ \\{zl By^W \dz\ 

(A18) 



We need to bound | log z|, ||(z/-A)-i||, and ||(z7-B)-i || 
on the outer and inner arcs 7j^ and 73 , and the upper and 
lower line segments 72 and 74 . 



We have | log z\ — | log \z\ + i arg(z) | < \/log^{R) + tt^ 
on 7*. On "fl we have 



|(z/-A)-i|| = |z-i|||(/-z-M)-i| 

00 

= \R-'\\\J2{z~'An 



k=0 



< |i?-i|£ii(z-M)ii 



(A19) 
(A20) 

(A21) 



fc=0 



< \z''\\{l~\\z-'A\\y' (A22) 



(i?-l)-i 



(A23) 



We used the fact that the above power series converges 
since ||z"^A|| = 1/R < 1. 

We cannot use the same argument on the 73 because 
||z~^A|| = R and the power series would diverge. There- 
fore, we rewrite {zl — A)^^ as 

(zI - A)-i = {zAA-^ - A)-i = -A-\I - zA-^y^ . 

(A24) 

The norm of A^^ is 1, so it suffices to show that the 
norm of (/ — zA^^)^^ bounded. To do this, we can use 
the power series expansion since ||zA^"'^|| = 1/R on 73. 
We obtain 



11(7- zA- 



< 



R 



R-1 



(A25) 



On 72 and 7I we have \\{zl — A)~^|| < 2 since ±e*5 
are the closest points to the imaginary axis, which could 
be in cr(A). The distance e'a^ and the imaginary axis is 
Ree*5 = 1/2. The same bounds are true for {zl — B). 

The length I7I of 7 is 7r(i?+l/i?) + 2(i?-l/i?^). Finally, 
we obtain the desired Lipschitz constant K by choosing 
it! = 2 which yields K < Gtt + 12. 



APPENDIX B: LIPSCHITZ CONTINUITY OF 
PARTITION FUNCTIONS WITH RESPECT TO 
HAMILTONIANS 

Let H and H be two self-adjoint operators with H = 
H + E and \\E\\ < f. Denote by Z{(3) and Z(/3) the cor- 
responding partition functions at (inverse) temperature 
p. We prove that 



{l/S)Z{P) < Z{I3) < SZ{f3) 



(Bl) 



holds for 6 — e'^'^/*. This makes it possible to estimate 
Z(P) with small relative error by estimating Z((3) instead 
provided that the norm of the error term E is sufficiently 
small for given temperature (3. 

This readily follows from Weyl's perturbation theorem 
[2]. Let A and B be two arbitrary self-adjoint operators 



7 



with eigenvalues \i{A) > . . . > A,\(A) and X[{B) > . . . > 
A^(i3), respectively. Then, 



max|A^^(yl)-Aj(B)| <\\A~B\ 



This implies 



(B2) 



(B3) 



We obtain the desired result by setting A = —j3H and 
B = -I3H. 



APPENDIX C: LIPSCHITZ CONTINUITY OF 
GIBBS STATES WITH RESPECT TO 
HAMILTONIANS 



Let H and H be two self-adjoint operators with H = 
H + E and \\E\\ < | where t = 7r/(4£',nax)- Denote by 
p{(3) = exp{-PH)/Z{(3) and = exp{-(3H)/Z{(3) 
the corresponding Gibbs states at (inverse) temperature 
p. We prove that the fidelity between p(/3) and p{(3) is 
high. More precisely, we have 

F{pi(3),pm>e~^'/'. (CI) 

Recall that the quantum fidelity is defined by 



(C2) 



F{P,P) = Tr J ^py/p. 



We bound the fildelity F{p, p) from below for the 
case of Gibbs states p = exp(— /3iJ)/Z(/3) and p = 
exp{—(3H)/Z{f3) by relying upon the following inequal- 
ity. It was proved in [S] that the function p i~> 



Tr(e^^/^eP^e^^/^) is increasing for self-adjoint opera- 
tors A and B inp e (0, oo). Its limit at p = is Tre^+'^. 
In particular, the inequalily 



Tre^+^ < Tr(ef^/2e^'^eP^/2 



(C3) 



holds for every p > 0. This is a strengthened variant of 
the Golden-Thompson inequality, which is obtained for 
p^l. 

Setting p ^ 2, A ^ ~f3H/2, and B = -~(3H/2, we 
bound the fidelity from below in terms of partition func- 
tions as follows 



Fip,p) 



rpj./'g-/3///2g-/3//g-/3ff/2 



1/2 



/Z(/3)Z(/3) 

rjj.(g-/3(H+H)/2> 



> 



Z(/3)Z(/3) 



(C4) 



(C5) 



(C6) 



/Z(/3)Z(/3) 

where H = {H + H]/2. Note that \\H - H\\ < e/{2t). 
Using the results of Appendix B, we obtain 



Fip,p)> 



Z(/3)i(/3) 



(C7) 



